Take Home Exercise 1

Published

January 30, 2023

Modified

February 12, 2023

In this exercise, we want to:

Conduct spatial point pattern analysis to discover the geographical distribution of functional and non-functional water points and their co-locations in Nigeria.

Setup

Loading Packages

pacman::p_load(sf, funModeling, maptools, raster, spatstat, tmap, tidyverse, sfdep)

Importing Data

Geospatial Data

NGA <- st_read("data/geospatial", layer = "nga_admbnda_adm2_osgof_20190417")%>%
  st_transform(crs = 26392)

Reading layer `nga_admbnda_adm2_osgof_20190417' from data source 
  `C:\zoe-chia\IS415\Take-home_Ex\Take-home_Ex01\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 774 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2.668534 ymin: 4.273007 xmax: 14.67882 ymax: 13.89442
Geodetic CRS:  WGS 84

Aspatial Data

Waterpoint Dataset

osun_wp_nga <- read_csv("data/aspatial/WPdx.csv") %>%
  filter(`#clean_country_name` == "Nigeria", 
         `#clean_adm1` == "Osun")

Data Handling

We want to combine our waterpoint aspatial data with our geospatial data. First, we have to give aspatial data its spatial properties. Then, we will clean our data before combining them with intersects().

Aspatial waterpoint data

Converting water point data into sf point features

Aspatial data does not have spatial properties. Hence, to turn it into geospatial data, we need to (i) convert it into the appropriate data type and (ii) add projection information (iii) transform to the respective coordinate system:

Convert the wkt field into sfc field by using st_as_sfc() data type.

osun_wp_nga$Geometry = st_as_sfc(osun_wp_nga$`New Georeferenced Column`)
osun_wp_nga

# A tibble: 5,745 × 75
   row_id `#source`      #lat_…¹ #lon_…² #repo…³ #stat…⁴ #wate…⁵ #wate…⁶ #wate…⁷
    <dbl> <chr>            <dbl>   <dbl> <chr>   <chr>   <chr>   <chr>   <chr>  
 1 225950 Federal Minis…    7.43    4.26 05/05/… Yes     Boreho… Well    Hand P…
 2 225524 Federal Minis…    7.78    4.56 04/22/… Yes     Protec… Well    Hand P…
 3 197014 Federal Minis…    7.49    4.53 04/30/… Yes     Boreho… Well    Mechan…
 4 225173 Federal Minis…    7.93    4.73 05/02/… Yes     Boreho… Well    Hand P…
 5 225843 Federal Minis…    7.74    4.44 05/08/… Yes     Boreho… Well    Hand P…
 6 235508 Federal Minis…    7.15    4.64 04/27/… Yes     Protec… Well    Hand P…
 7 197708 Federal Minis…    7.87    4.72 05/13/… Yes     Boreho… Well    Mechan…
 8 195041 Federal Minis…    7.73    4.45 06/17/… Yes     Protec… Spring  <NA>   
 9 225222 Federal Minis…    7.81    4.15 05/14/… Yes     Protec… Spring  Mechan…
10 460770 GRID3             7.4     4.33 06/13/… Unknown Boreho… Well    <NA>   
# … with 5,735 more rows, 66 more variables: `#water_tech_category` <chr>,
#   `#facility_type` <chr>, `#clean_country_name` <chr>, `#clean_adm1` <chr>,
#   `#clean_adm2` <chr>, `#clean_adm3` <chr>, `#clean_adm4` <chr>,
#   `#install_year` <dbl>, `#installer` <chr>, `#rehab_year` <lgl>,
#   `#rehabilitator` <lgl>, `#management_clean` <chr>, `#status_clean` <chr>,
#   `#pay` <chr>, `#fecal_coliform_presence` <chr>,
#   `#fecal_coliform_value` <dbl>, `#subjective_quality` <chr>, …

Convert tibble dataframe into an sf object using st_sf(). It is also important for us to include the referencing system of the data into the sf object. Data in wp_nga is in decimals, and in tibble form. It does not have the projection information. We need to put back its original projection inform (WGS84).

osun_wp_sf <- st_sf(osun_wp_nga, crs = 4326)
osun_wp_sf

Simple feature collection with 5745 features and 74 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 4.032004 ymin: 7.060309 xmax: 5.06 ymax: 8.061898
Geodetic CRS:  WGS 84
# A tibble: 5,745 × 75
   row_id `#source`      #lat_…¹ #lon_…² #repo…³ #stat…⁴ #wate…⁵ #wate…⁶ #wate…⁷
 *  <dbl> <chr>            <dbl>   <dbl> <chr>   <chr>   <chr>   <chr>   <chr>  
 1 225950 Federal Minis…    7.43    4.26 05/05/… Yes     Boreho… Well    Hand P…
 2 225524 Federal Minis…    7.78    4.56 04/22/… Yes     Protec… Well    Hand P…
 3 197014 Federal Minis…    7.49    4.53 04/30/… Yes     Boreho… Well    Mechan…
 4 225173 Federal Minis…    7.93    4.73 05/02/… Yes     Boreho… Well    Hand P…
 5 225843 Federal Minis…    7.74    4.44 05/08/… Yes     Boreho… Well    Hand P…
 6 235508 Federal Minis…    7.15    4.64 04/27/… Yes     Protec… Well    Hand P…
 7 197708 Federal Minis…    7.87    4.72 05/13/… Yes     Boreho… Well    Mechan…
 8 195041 Federal Minis…    7.73    4.45 06/17/… Yes     Protec… Spring  <NA>   
 9 225222 Federal Minis…    7.81    4.15 05/14/… Yes     Protec… Spring  Mechan…
10 460770 GRID3             7.4     4.33 06/13/… Unknown Boreho… Well    <NA>   
# … with 5,735 more rows, 66 more variables: `#water_tech_category` <chr>,
#   `#facility_type` <chr>, `#clean_country_name` <chr>, `#clean_adm1` <chr>,
#   `#clean_adm2` <chr>, `#clean_adm3` <chr>, `#clean_adm4` <chr>,
#   `#install_year` <dbl>, `#installer` <chr>, `#rehab_year` <lgl>,
#   `#rehabilitator` <lgl>, `#management_clean` <chr>, `#status_clean` <chr>,
#   `#pay` <chr>, `#fecal_coliform_presence` <chr>,
#   `#fecal_coliform_value` <dbl>, `#subjective_quality` <chr>, …

Transform to Nigeria’s projected coordinate system

osun_wp_sf <- osun_wp_sf %>%
  st_transform(crs = 26392)

Data Wrangling for Water Point Data

Get status_clean column

Rename ‘#status_clean’ column and replace NA with ‘unknown’.

osun_wp_sf_nga <- osun_wp_sf %>%
  dplyr::rename(status_clean = '#status_clean') %>%
  dplyr::select(status_clean) %>%
  dplyr::mutate(status_clean = replace_na(
    status_clean, "unknown"
  ))

Visualise status_clean column

funModeling::freq(data = osun_wp_sf_nga, 
     input = 'status_clean')

               status_clean frequency percentage cumulative_perc
1                Functional      2406      41.88           41.88
2            Non-Functional      2086      36.31           78.19
3                   unknown       748      13.02           91.21
4  Functional, needs repair       259       4.51           95.72
5       Non-Functional, dry       159       2.77           98.49
6    Functional, not in use        64       1.11           99.60
7  Abandoned/Decommissioned        15       0.26           99.86
8 Functional but not in use         8       0.14          100.00

Group data by filter() to obtain functional, non-functional and unknown statuses

wp_functional <- osun_wp_sf_nga %>% 
  filter(status_clean %in% 
           c("Functional",
             "Functional but not in use", 
             "Functional but needs repair", 
             "Functional, needs repair"))


wp_nonfunctional <- osun_wp_sf_nga %>% 
  filter(status_clean %in%
           c("Abandoned/Decommissioned",
             "Abandoned",
             "Non-Functional due to dry season", 
             "Non-Functional",
             "Non functional due to dry season", 
             "Non-Functional, dry"))

wp_unknown <- osun_wp_sf_nga %>% 
  filter(status_clean == "unknown")

funModeling::freq(data = wp_functional, 
     input = 'status_clean')

               status_clean frequency percentage cumulative_perc
1                Functional      2406      90.01           90.01
2  Functional, needs repair       259       9.69           99.70
3 Functional but not in use         8       0.30          100.00

Add a new column to categorise water points’ functionality

osun_wp_sf_nga$new_status <- NA
osun_wp_sf_nga$new_status[osun_wp_sf_nga$status_clean %in% c("Functional",
             "Functional, not in use", 
             "Functional but not in use",
             "Functional but needs repair",
             "Functional, needs repair")] <- "Functional"
osun_wp_sf_nga$new_status[osun_wp_sf_nga$status_clean %in%
           c("Abandoned/Decommissioned",
             "Abandoned",
             "Non-Functional due to dry season", 
             "Non-Functional",
             "Non functional due to dry season", 
             "Non-Functional, dry")] <- "Non-Functional"
osun_wp_sf_nga$new_status[osun_wp_sf_nga$status_clean %in% c("unknown", "NA")] <- "Unknown"

funModeling::freq(data = osun_wp_sf_nga, 
     input = 'new_status')

      new_status frequency percentage cumulative_perc
1     Functional      2737      47.64           47.64
2 Non-Functional      2260      39.34           86.98
3        Unknown       748      13.02          100.00

Geospatial Data

Get Study Area, Osun State

NGA <- NGA %>% 
  filter(`ADM1_EN` == "Osun")

Check for duplicates

NGA$ADM2_EN[duplicated(NGA$ADM2_EN)==TRUE]

character(0)

There are no duplicates.

Point in Polygon Count

To find out the number of total, functional, nonfunctional and unknown water points.

NGA_wp <- NGA %>% 
  mutate(`total_wp` = lengths(
    st_intersects(NGA, osun_wp_sf_nga))) %>%
  mutate(`wp_functional` = lengths(
    st_intersects(NGA, wp_functional))) %>%
  mutate(`wp_nonfunctional` = lengths(
    st_intersects(NGA, wp_nonfunctional))) %>%
  mutate(`wp_unknown` = lengths(
    st_intersects(NGA, wp_unknown)))

Visualising

ggplot(data = NGA_wp,
       aes(x = total_wp)) + 
  geom_histogram(bins=20,
                 color="black",
                 fill="light blue") +
  geom_vline(aes(xintercept=mean(
    total_wp, na.rm=T)),
             color="red", 
             linetype="dashed", 
             size=0.8) +
  ggtitle("Distribution of total water points by LGA") +
  xlab("No. of water points") +
  ylab("No. of\nLGAs") +
  theme(axis.title.y=element_text(angle = 0))

Save data in rds format

write_rds(NGA_wp, "data/rds/NGA_wp.rds")

Import the saved data in rds

NGA_wp <- read_rds("data/rds/NGA_wp.rds")

Exploratory Data Analysis

Visualising Distribution of Water Points with Choropleth Maps

p1 <- tm_basemap(server="OpenStreetMap") + 
  tm_shape(NGA_wp) +
  tm_fill("wp_functional",
          n = 10,
          style = "equal",
          palette = "Blues") +
  tm_borders(lwd = 0.1,
             alpha = 1) +
  tm_layout(main.title = "Distribution of functional water point by cities",
            legend.outside = FALSE,
            main.title.size = 1) + 
  tm_view(set.zoom.limits = c(8, 10))

p2 <- tm_basemap(server="OpenStreetMap") +
  tm_shape(NGA_wp) +
  tm_fill("wp_nonfunctional",
          n = 10,
          style = "equal",
          palette = "Blues") +
  tm_borders(lwd = 0.1,
             alpha = 1) +
  tm_layout(main.title = "Distribution of non-functional water point by cities",
            legend.outside = FALSE,
            main.title.size = 1) + 
  tm_view(set.zoom.limits = c(8, 10))

tmap_mode('view')
tmap_arrange(p1, p2, nrow = 1, widths = c(0.5, 0.5))

tmap_mode('plot')

tmap_mode('view')
tm_basemap(server="OpenStreetMap")+
tm_shape(NGA) + 
  tm_polygons() + 
  tm_shape(wp_functional) + 
  tm_dots(size = 0.01, 
          col = "yellow", 
          border.lwd = 0.5, 
          legend.show = TRUE
          ) + 
  tm_layout(legend.outside = FALSE) + 
  tm_shape(wp_nonfunctional) + 
  tm_dots(col="turquoise",
    size = 0.01, 
          border.col = "red", 
          border.lwd = 0.5,
    legend.show = TRUE) +
  tm_layout(legend.outside = FALSE)+
  tm_view(set.zoom.limits = c(9, 12))

tmap_mode('plot')

Converting sf data frames to sp’s Spatial* class

osun_wp_spdf <- as_Spatial(osun_wp_sf) 
osun_wp_nga_spdf <- as_Spatial(NGA_wp)

osun_wp_spdf

class       : SpatialPointsDataFrame 
features    : 5745 
extent      : 177285.9, 291287.1, 340054.1, 450859.7  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
variables   : 74
names       : row_id,                                     X.source, X.lat_deg, X.lon_deg,          X.report_date, X.status_id,    X.water_source_clean, X.water_source_category, X.water_tech_clean, X.water_tech_category, X.facility_type, X.clean_country_name, X.clean_adm1, X.clean_adm2, X.clean_adm3, ... 
min values  : 154218, Federal Ministry of Water Resources, Nigeria,  7.060309, 4.0320038, 01/01/2010 12:00:00 AM,          No,                Borehole,                  Spring,          Hand Pump,             Hand Pump,        Improved,              Nigeria,         Osun,     Aiyedade,           NA, ... 
max values  : 685450,                                        GRID3, 8.0618983,      5.06, 09/16/2015 12:00:00 AM,         Yes, Undefined Hand Dug Well,                    Well,           Tapstand,              Tapstand,        Improved,              Nigeria,         Osun,       Osogbo,           NA, ...

osun_wp_nga_spdf

class       : SpatialPolygonsDataFrame 
features    : 30 
extent      : 176503.2, 291043.8, 331434.7, 454520.1  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
variables   : 20
names       :    Shape_Leng,       Shape_Area,  ADM2_EN, ADM2_PCODE, ADM2_REF, ADM2ALT1EN, ADM2ALT2EN, ADM1_EN, ADM1_PCODE, ADM0_EN, ADM0_PCODE,  date, validOn, validTo,        SD_EN, ... 
min values  : 0.26445678806, 0.00248649736648, Aiyedade,   NG030001, Aiyedade,         NA,         NA,    Osun,      NG030, Nigeria,         NG, 17134,   18003,      NA, Osun Central, ... 
max values  :  1.8470166597,  0.0737271661922,   Osogbo,   NG030030,   Osogbo,         NA,         NA,    Osun,      NG030, Nigeria,         NG, 17134,   18003,      NA,    Osun West, ...

Converting the Spatial* class into generic sp format

osun_wp_sp <- as(osun_wp_spdf, "SpatialPoints")
osun_wp_nga_sp <- as(osun_wp_nga_spdf, "SpatialPolygons")

Converting to PPP format

osun_wp_ppp <- as(osun_wp_sp, "ppp")

Converting non-functional and functional water points from sf format to ppp format

osun_functional_spdf <- as_Spatial(wp_functional)
osun_functional_sp <- as(osun_functional_spdf, "SpatialPoints")
osun_functional_ppp <- as(osun_functional_sp, "ppp")

osun_nonfunctional_spdf <- as_Spatial(wp_nonfunctional)
osun_nonfunctional_sp <- as(osun_nonfunctional_spdf, "SpatialPoints")
osun_nonfunctional_ppp <- as(osun_nonfunctional_sp, "ppp")

par(mfrow=c(1,3))
plot(osun_wp_ppp, main = "Osun Waterpoints (PPP)")
plot(osun_functional_ppp, main = "Functional Waterpoints (PPP)")
plot(osun_nonfunctional_ppp, main = "Non-Functional Waterpoints (PPP)")

Check for duplicated points

any(duplicated(osun_wp_ppp))

[1] FALSE

Create owin object

osun_owin <- as(osun_wp_nga_sp, "owin")
plot(osun_owin, main = "Osun")

Combining point events object and owin object

osun_ppp = osun_wp_ppp[osun_owin]
functional_osun_owin_ppp = osun_functional_ppp[osun_owin]
nonfunctional_osun_owin_ppp = osun_nonfunctional_ppp[osun_owin]

Transforming unit measurements from metre to kilometre

osun_ppp.km = rescale(osun_ppp, 1000, "km")
functional_osun_owin_ppp.km = rescale(functional_osun_owin_ppp, 1000, "km")
nonfunctional_osun_owin_ppp.km = rescale(nonfunctional_osun_owin_ppp, 1000, "km")

par(mfrow=c(1,3))
plot(osun_ppp, main = "Osun PPP")
plot(functional_osun_owin_ppp, main = "Functional WP PPP")
plot(nonfunctional_osun_owin_ppp, main = "Non-Functional WP PPP")

Kernel Density Estimation (KDE)

Computing KDE of waterpoints in Osun using automatic and fixed bandwidth selection method

Adaptive bandwidth

kde_osunwp_abw <- density(osun_ppp.km,
                         sigma = bw.diggle,
                         edge = TRUE, 
                         kernel = "gaussian")

kde_functional_abw <- density(functional_osun_owin_ppp.km,
                              sigma = bw.diggle,
                              edge = TRUE, 
                              kernel = "gaussian")

kde_nonfunctional_abw <- density(nonfunctional_osun_owin_ppp.km,
                              sigma = bw.diggle,
                              edge = TRUE, 
                              kernel = "gaussian")

kde_osunwp_abw

real-valued pixel image
128 x 128 pixel array (ny, nx)
enclosing rectangle: [176.5, 291.04] x [331.43, 454.52] km

Fixed Bandwidth

kde_osunwp_fbw <- density(osun_ppp.km,
             sigma = 0.6, 
             edge = TRUE, 
             kernel = "gaussian")

kde_functional_fbw <- density(functional_osun_owin_ppp.km,
                              sigma = 0.6,
                              edge = TRUE, 
                              kernel = "gaussian")

kde_nonfunctional_fbw <- density(functional_osun_owin_ppp.km,
                              sigma = 0.6,
                              edge = TRUE, 
                              kernel = "gaussian")

Comparing Fixed and Adaptive KDEs:

par(mfrow=c(3,2))
plot(kde_osunwp_abw, main = "Total WP KDE Adaptive Bandwidth")
plot(kde_osunwp_fbw, main = "Total WP KDE Fixed Bandwidth")
plot(kde_functional_abw, main = "Functional WP KDE Adaptive Bandwidth")
plot(kde_functional_fbw, main = "Functional WP KDE Fixed Bandwidth")
plot(kde_nonfunctional_abw, main = "Non-Functional WP KDE Adaptive Bandwidth")
plot(kde_nonfunctional_fbw, main = "Non-Functional WP KDE Fixed Bandwidth")

Converting KDE Output into grid and raster object

For mapping purposes:

# Total WPs 
gridded_kde_osunwp_fbw <- as.SpatialGridDataFrame.im(kde_osunwp_fbw)
kde_osunwp_fbw_raster <- raster(gridded_kde_osunwp_fbw)

# Functional WPs
gridded_kde_functional_fbw <- as.SpatialGridDataFrame.im(kde_functional_fbw)
kde_functional_fbw_raster <- raster(gridded_kde_functional_fbw)

#Non-Functional WPs
gridded_kde_nonfunctional_fbw <- as.SpatialGridDataFrame.im(kde_nonfunctional_fbw)
kde_nonfunctional_fbw_raster <- raster(gridded_kde_nonfunctional_fbw)

Plot KDE in OpenStreeMap

tmap_mode("view")
# Total WP
totalwpraster <- tm_shape(kde_osunwp_fbw_raster) + 
  tm_basemap(server="OpenStreetMap") + 
  tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE, title = "Total WP Raster")+
  tm_view(set.zoom.limits = c(18,20))

# Functional WP
functional_raster <- tm_shape(kde_functional_fbw_raster) + 
  tm_basemap(server="OpenStreetMap") + 
  tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE, title = "Functional WP Raster")+
  tm_view(set.zoom.limits = c(18,20))

# Non-Functional WP
nonfunctional_raster <- tm_shape(kde_nonfunctional_fbw_raster) +
  tm_basemap(server="OpenStreetMap") + 
  tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE, title = "Non-Functional WP Raster")+
  tm_view(set.zoom.limits = c(18,20))

tmap_arrange(totalwpraster, functional_raster, nonfunctional_raster, nrow = 1)

tmap_mode('plot')

Conclusions

Spatial Pattern Observations

There seems to be very slight clustering for the functional and non-functional water points. From the Kernel Density maps, they seem to be clustered around the same areas. However, overall, they seem to be randomly distributed.

Advantage of Kernel Density Map over Point Map

Kernel Density map smooths out noise and highlights any clusters that may be present. It shows us the separation between feature and noise points. However, the point map shows us all points including noise which may be difficult to identify clear clusters.

Nearest Neighbour Analysis

The test hypotheses are:

Ho = The distribution of water points are randomly distributed.

H1 = The distribution of water points are not randomly distributed.

The 95% confident interval will be used.

Clark and Evans Test

clarkevans.test(osun_ppp,
                correction="none",
                clipregion="osun_owin", 
                alternative=c("two.sided"),
                nsim=99)


    Clark-Evans test
    No edge correction
    Monte Carlo test based on 99 simulations of CSR with fixed n

data:  osun_ppp
R = 0.40209, p-value = 0.02
alternative hypothesis: two-sided

Conclusions

With the p-value of 0.02 being less than the critical level of 0.05, we should reject the null hypothesis and conclude that the distribution of water points is not random.

Analysing Spatial Point Process Using L-Function

All water points

L Function Estimation

#L_wp = Lest(osun_ppp, correction = "Ripley")

#plot(L_wp, . -r ~ r, ylab= "L(d)-r", xlab = "d(m)")

Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of water points at Osun State are randomly distributed.

H1 = The distribution of water points at Osun State are not randomly distributed.

The null hypothesis will be rejected if p-value if smaller than alpha value of 0.05.

The code chunk below is used to perform the hypothesis testing at the 95% confidence level.

#L_wp.csr <- envelope(osun_ppp, Lest, nsim = 39, rank = 1, glocal=TRUE)

# plot(L_wp.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

Functional WP

L Function Estimation

# L_ck = Lest(functional_osun_owin_ppp, correction = "Ripley")

# plot(L_ck, . -r ~ r, ylab= "L(d)-r", xlab = "d(m)")

Complete Spatial Randomness Test

# L_ck.csr <- envelope(functional_osun_owin_ppp, Lest, nsim = 39, rank = 1, glocal=TRUE)

# plot(L_ck.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

Non-Functional WP

L Function Estimation

# L_nonfwp = Lest(nonfunctional_osun_owin_ppp, correction = "Ripley")

# plot(L_nonfwp, . -r ~ r, ylab= "L(d)-r", xlab = "d(m)")

Complete Spatial Randomness Test

# L_nonfwp.csr <- envelope(functional_osun_owin_ppp, Lest, nsim = 39, rank = 1, glocal=TRUE)

# plot(L_nonfwp.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

Unable to run L Function

I was unable to run the L function in time, but I will update this section with another desktop and input the results here.

Local Colocation Quotients

Visualising the sf layers

tmap_mode("view")
tm_shape(NGA_wp) +
  tm_polygons() +
tm_shape(osun_wp_sf_nga)+ 
  tm_dots(col = "new_status",
             size = 0.01,
             border.col = "black",
             border.lwd = 0.5) +
  tm_view(set.zoom.limits = c(9, 10))

tmap_mode("plot")

Preparing nearest neighbours list

nb <- include_self(
  st_knn(st_geometry(osun_wp_sf_nga), 6))

Compute kernel weights

wt <- st_kernel_weights(nb, osun_wp_sf_nga, "gaussian", adaptive= TRUE)

Preparing the vector list

Functional <- osun_wp_sf_nga %>%
  filter(new_status == "Functional")
target <- Functional$new_status

NonFunctional <- osun_wp_sf_nga %>%
  filter(new_status == "Non-Functional")
B <- NonFunctional$new_status

Run 50 simulations:

LCLQ <- local_colocation (target, B, nb, wt, 49)

Combine waterpoint and local location table

LCLQ_wp <-cbind(osun_wp_sf_nga, LCLQ)

tmap_mode('view')
tm_shape(NGA_wp) + 
tm_polygons() + 
  tm_shape(LCLQ_wp) + 
  tm_dots(col = "Non.Functional", 
          size = 0.01, 
          border.col = "black", 
          border.lwd = 0.5) + 
  tm_view(set.zoom.limits = c(9, 12))

tmap_mode('plot')

Conclusions

From the results above, we can see that there are many missing, or undefined LCLQ types. This means that the functional water points do not have other features within its neighbourhood, and hence do not have any correlation and are independent from each other.